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^^ Abstract In this work, a geometric discretization of the Navier-Stokes equations 

lO is sought by treating momentum as a covector-valued volume-form. The novelty 

^N of this approach is that we treat conservation of momentum as a tensor equation 

i__i and describe a higher order approximation to this tensor equation. The resulting 

.^ scheme satisfies mass and momentum conservation laws exactly, and resembles a 

staggered-mesh finite-volume method. Numerical test-cases to which the discretiza- 
tion scheme is applied are the Kovasznay flow, and lid-driven cavity flow. 






X 



1 Navier-Stokes equations 



K^ Mimetic discretizations aim to represent physics in a discrete sense, in contrast to 

,__! differential formulations, which are concerned with the limit h — > 0. For the case in 

(^\ which /i 7^ geometrical considerations play an important role in the correct discrete 

'3\ formulation, m [U |4l |6l [TOl . Application of these ideas to continuum models are 

^^ described in fSl, Appendix A] and |8 13 1. The novel aspect in this paper is that 

t:J" continuum ideas are applied to incompressible, viscous flows using spectral basis 

^D functions. 

^^ We start with the incompressible Navier-Stokes equations (p — 1), written in the 

. . integral formulation, as given in many textbooks and we try to make precise what 

these statements mean. It is important to give an accurate meaning to all variables, 

because when we want to represent these physical quantities on finite grids, we 

'^ want to preserve the main structure of the equations. Conservation of mass (p — 1) 
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is usually given by 



vnd5 = 0, (1) 

dn 



and conservation of momentum. 



I 

Jf. 



-, dV+ / \(E)\-ndS= / andS (2) 

n at Jdn Jdn 

and Newtonian stress relation 

CJ = -pI + Ai(Vv+(Vv)^) . (3) 

Here \,p,CT and jj. denote velocity, pressure, total stress tensor and dynamic vis- 
cosity, respectively; I and n are the identity matrix and the outward unit normal to 
the boundary, respectively. The above are balance equations for volumetric quanti- 
ties that depend on their fluxes through surfaces and are more physical than their 
differential counterparts. 



1.1 Momentum and velocity 

The first term in (J2]) indicates that velocity (and its time derivative) can be integrated 
over a volume. But velocity is generally not associated to volumes, but is defined 
as the tangent vector at a given point along the trajectory of a particle. Velocity is 
therefore a vector-valued 0-form. This statement means that to every point in space- 
time (a zero-dimensional object) we associate a vector Let V be the linear vector- 
space of all possible vectors at a given point in space, then we can define the space 
V* of all linear functionals on V. Elements of V* are called covectors. The spaces 
V and V* are isomorphic, but there is no canonical isomorphism which relates an 
element v e V to an element a G V*. Once a metric is defined, one can associate 
with every vector at a point a corresponding covector. This map is called the flat 
operator: b : V — > V*. The covector associated with a vector v is then denoted by v" . 

The linear vector space V associated to a point p is called the tangent space at p, 
denoted by TpQ . The corresponding dual space is caUed the cotangent space at p 
denoted by T*£l. The collection of all tangent spaces in the domain CI is called the 
tangent bundle, TQ and the collection of cotangent spaces is called the cotangent 
bundle, T*Q. Let a G T*Q and v G TQ, then (o:,v) associates to each point p in 
Q the value a\p{\\p). 

With every ^-form we can associate a (n — A:)-form with a different type of orien- 
tation, see L2J. The collection of all fc-forms on Q is denoted by A*^(i2). The metric 
dependent operator which establishes this connection is the Hodge-* operator For 
continuum models we need to combine the b and Hodge-* into the operator •^ (see 
also ifTSJI for such operations) 



*^ : TQ ® A*(f2) -^ T*Q (g)A"-''{Q) 
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If we apply this operator to velocity v G TQ (>^A^{Q) we obtain 

Similarly, we can define •f : T*Q <SiA''{Q) -^ TQ <SiA"-''{Q). The physical quan- 
tity m is called momentum density or the momentum per unit volume. This is a 
covector-valued volume form. So instead of integrating 'velocity' over the domain 
we are tempted to write 

Jn Jn 
This integral is not defined, because it assumes that we can integrate over the tangent 
spaces in Q. The basis in each tangent space, however, may differ from point to 

point. In order to define the momentum integral we introduce the operator A 

A: (T*Q(x)A''{Q)]®(Tn®A'{Q)] ^A*+'(i2), 
given by a e T*Q<»A''{Q) and we TQ<g)A'{Q) 

a A w = (a,w)djc'*^' Adjc*'' . 

This operation yields a (^+Z)-form which can be integrated over (A; +Z) -dimensional 
submanifolds. 

If we apply momentum density m to any vector field w (not necessarily a velocity 

field) using this operator we get m A w G A"(i2) and this can be integrated over a 
volume. So the proper way to interpret the time rate of change of momentum should 
be 

/ ^•^(v)Aw, Vweri2(8)A*'(i2). (4) 

Jn at 

In many textbooks on fluid the distinction between dynamics momentum density 
(usually called 'momentum') and velocity is ignored; one is just a scalar multiple 
of the other, m =^ pv, but the use of the vector w in (HI is generally incorporated. 
The textbooks then say: 'We consider this equation for each component separately 
..'. This is a strange sentence, because components have no physical relevance, only 
vectors, i.e. components plus associated basis vectors are physically relevant. But 
what is meant by this statement is that for the vector field w in Q a uniform vector 
field in the jc'-direction is taken. The generality 'all vector fields' is in these text- 
books compensated by the fact that momentum conservation should hold for 'all 
volumes'. 
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1.2 Convection 

Now that we understand how momentum density should be integrated over a vol- 
ume, we can also define convection of momentum density. After pairing with an 
arbitrary vector field, w, we obtain a volume form and we apply the Lie derivative 
to this volume form, see IHTII . The Lie derivative for a volume form, j3^"', is given 
by 

=2'v/3W=divj3W, 

and then the generalized Stokes theorem converts this exact form to a boundary 
integral 



^v'«Aw= / iv(mAw). (5) 

n Jdii 

Compare this expression with the the convective term in (J2]) and note that it does 
not require an inner product nor the definition of an outward unit normal. The in- 
ner product is avoided since we work with differential forms and duality pairing is 
metric-free and the orientation of the elements in the mesh, ||6l, avoids the use of 
explicitly defined normals. 



1.3 Stress tensor and surface force density 

The last term in (|2]i denotes the action of the viscous forces on the flow represented 
by the stress tensor a. The stress tensor is an infinitesimal quantity in the limit 
for h -^ 0. On a finite mesh we can identify volumes over which we integrate the 
momentum density and the boundary of these volumes where surface forces act. 
In continuum mechanics forces are 'smeared out', so we introduce the surface force 
density givenhy tG T*Q^A"^^{Q). This is acovector-valued {n~ \)-form. Forces 
are generally associated with covectors, (2, JL2|, and in the current setting need to 
be covectors in order to equate them to the time rate of change of momentum which 
was also covector-valued. It is furthermore a («— l)-form since it acts on the the 
boundary of n-dimensional volumes, see also fs! Appendix A] and ll8l [T3]| . Again, 
covector-valued forms cannot be integrated, so the proper way is pair this to pair 
it with an arbitrary vector field w before integration over surfaces is possible. The 
momentum equation then becomes 



d 
dt 



/•^(v)Aw+/ iv(*^(v)Aw)= / tAw, VwerX2®A''(X2). (6) 
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1.4 Newtonian stress relation 

The pressure scalar is an outer-oriented volume form, p'"'. Pressure force density is 
represented as a covector- valued (n — l)-form 

p== (*/?) djc'^djc' A...djt:'...Adjt:" , 

where the notation^indicates that this term is omitted and dx' di' A . . . dx' . . . A di" 
is the identity tensor, see also example |J5] §9. 3a]. This description agrees with ||9] 

for Stokes flow. Note that p A w = iwp'"'- 

The velocity gradient is represented as the covariant differential of the velocity 
vector field, Vv which is a vector-valued 1-form, see |5, §9. 3b]. In this paper we 
restrict ourselves to Euclidean space for which the connection 1 -forms vanish. Ap- 
plying *u(Vv) transforms the vector- valued 1-form into a covector-valued (n — 1)- 
form, where the diffusion coefficient is contained in the Hodge-* operator. In this 
paper we assume /i to be constant. 



1.5 Conservation of mass 

Let co'"' be the standard volume form, then the divergence of a vector field is defined 
as (divv)a)'"' = .ifyO^"' — dx^a^"' . Integration over a volume and applying Stokes 
theorem gives 

Jn Jdii 

This is the proper translation of ([T} as found in textbooks on incompressible flow. 
The velocity flux field, ivfi)'"\ is isomorphic to the velocity vector field. The velocity 
flux field will be used in the discrete representation of the Navier-Stokes equations. 
The relation between the velocity fluxes iyfi)'"' and * {\) is given by 

•^ (v) A w = w^ A ivOJ^"* , ^weTQ^A^iQ) . (!) 

The volume forms and (n — 1) -forms appearing in all integrals are all outer- 
oriented. 



2 Discrete representation 

In the full differential geometric setting as described above, the integration only 
makes sense when paired with all vector fields w. Here we choose the uniform 
vector field in the x- and y-direction only and impose that conservation should hold 
for all volumes in our spectral elements. These volumes are generated by the Gauss- 
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Lobatto grid in the spectral element and will be denoted by Qjj. So in this section w 
is either d^ or dy. 

Figure [T] displays one spectral element and its Gauss-Lobatto grid (solid lines) in 
2D. The dotted gray Unes represent the dual grid, see Ii6,il0j. 
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Fig. 1 The velocities are discretized as outer-oriented mass-fluxes, and live on surfaces (S) of 
the Gauss-Lobatto grid shown above, while pressure is discretized on volumes (Q). Momenta 
are discretized on staggered volumes (Q) and their fluxes on the surfaces (5) surrounding these 
staggered volumes. 



Momentum is reduced onto a volume consisting of a primal (n — 1) -chain and 
a dual 1 -chain. In 2D these volumes consist of tensor products of primal and dual 
edge as shown in Figure [Tlby volumes enclosed by solid (primal) and dashed (dual) 
lines. The location of the unknowns coincides with those in staggered finite volume 
methods. The difference is that in this formulation the unknowns represent integral 
values, whereas in finite volume methods the unknowns either represent average or 
nodal values. Let us denote the primal surfaces by St, then discrete velocity is given 
by 



Si 

This yields a metric -free description of conservation as mass as shown in ifTOl lTJI. 
The reduction of the pressure field is on outer-oriented volumes, see also fTO^. 

Integrals of momentum flux, ^^"^^\ pressure force, iw/?'"', and velocity gradi- 
ents are represented on the boundary of the momentum volumes indicated in Fig- 
ure [U 
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Once we have the discrete variables for mass flux, momentum and pressure, we 
use the spectral element functions described in ifTOl lTl, to interpolate these values in 
such a way that the integral values are preserved. 

Using (|7]i we can write the relation between momentum and velocity flux as 

/ mAw- / w^Aiv©'"* =0 — >my,-P'^u = , (8) 

Jntj Jhij 

where tljj are the volume where momentum is reduced, see Fig. Ill and P™ is the 
matrix which maps discrete velocity (which is discretized as mass-fluxes) to discrete 
momentum (on the staggered-grid). The discrete representation of momentum-flux, 
pressure force iwp'"' and traction forces, ^^ (V„v) (which can be equivalently writ- 
ten as Tk-^ (V*" •"v) A w = •„ (Vt^H m) A w, which, in Cartesian coordinates and with 
constant w becomes d* {m A w)), are given by 

• Convective-flux, see [111, ,^^ = iy{m A w): 

f^i'^j3(')) -fmAw,v^Aj3(')) , = 0. 

\ Jn \ Jn (9) 

— >Mi 1 #w - Cyihv, = . 

• Pressure-force, J^ ~i-„p^^': 

J^y^ - p(2) (w) = ^^4 - P^P = Bp . (10) 

1) 



• Diffusive-fluxes, ^ ' =d*Jm f\y/) 



V Jn \ Jn Jdn 



>Mn^y,-DlM22m^BT 



(11) 



The discrete continuity equation is given by 

DziM^O. (12) 

In the above, j3^'' is an arbitrary 1-form; Mn and M22 are mass-matrices for 1- 
and 2-forms on the staggered mesh; Cv is the convection matrix and depends on v 
(which can be retrieved from reconstruction of v using the edge basis, |6|); P^ is 
the matrix which converts the scalar p to pressure-force 1-forms; Bp and Bp are 
the boundary integrals for pressure and stress, respectively, obtained from integra- 
tion by parts; and D21 and D21 are incidence matrices which discretely represent 
the exterior-derivative with entries containing only { — 1,0, 1}. The algebraic system 
thus obtained is solved for v and p for w = {dx,dy}. 
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(a) Pressure, /j-refinement 



(b) Pressure, p-refinement 




(c) Velocity, /i-refinement 



(d) Velocity, p-refinement 



Fig. 2 Convergence plots for Kovasznay flow with mesh size, h, and order, p. Optimal rates are 
shown (^^) for the /i-refinement cases. nElemX and nElemY refer to the number of elements in 
X and Y directions, and p refers to the order of elements used. 



3 Results 



3.1 Kovasznay Flow 



Kovasznay flow is an analytical solution to Navier-Stokes' equations. The so- 
lution is M = 1 — e cos{27ty), v — ^e sin(2ny) and p = ^(l — e ), where 



2v " 



1 

4v2 



A'K?-. The kinematic-viscosity chosen for this flow was V = ^ 



and the computational domain considered was Q. = [—0.5 1] x [—0.5 —0.5]. The 
/!,/7-adaptivity plots for this problem are given in Fig.l2]for pressures (Fig. |2(a)] and 
Fig. 2(b) I and velocities (Fig. 2(c) and Fig. 2(d) i. It can be seen that the solutions 



converge exponentially and optimally. There is some stagnation observed in con- 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



(a) Stream-function contours 



(b) Pressure contours 




(c) X-velocity at x = 0.5 



(d) Y-velocity at y = 0.5 



Fig. 3 (Top) Streamfunction and Pressure contours with a single spectral element of order 16. 
(Bottom) Centerline velocities are plotted (^^) and compared with the solutions of 1 3 1 (^^), and 
the solutions are found to be reasonably close. Mesh size is 4 x 4 and made up of elements of order 
6. 



vergence for a mesh with a single element, and this is attributed to the fact that our 
basis may not be capturing certain modes (even/odd). 



3.2 Lid-driven Cavity Flow 



The second numerical test-case chosen was the classic Lid-driven cavity flow on 
a unit square domain with the top-lid velocity, m/, = — 1 and a Reynolds number 
of 1000. The solutions for the pressure and streamfunction contours calculated 
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for a single spectral element of order p = 16 are shown in the top-half of Fig. [3] 
Centerline-velocity solutions with a lower order of p — 6 but with multiple ele- 
ments (4x4 mesh) and comparisons with the results of [3| are also shown in the 
bottom-half of Fig. [3] Good agreement is seen between the benchmark results and 
our results. 



References 

[1] P B Bochev and J M Hyman. Principles of mimetic discretizations of differ- 
ential operators. IMA Volumes In Mathematics and its Applications, 142:89, 

2006. 
[2] A. Bossavit. Handbook of Numerical Analysis, volume 13, chapter Discretiza- 
tion of electromagnetic problems, pages 105-197. Elsevier, 2005. 
[3] O. Botella and R. Peyret. Benchmark spectral results on the lid-driven cavity 

flow. Computers & Fluids, 27(4):421-433, 1998. 
[4] M Desbrun and AN Hirani. Discrete exterior calculus. Arxiv preprint 

arXiv:math/0508341 , 2005. URL |http: //arxiv. org/abs/math71 

105083411 
[5] Th. Frankel. The geometry of physics: An introduction. Cambridge University 

Press, 2011. 
[6] M Gerritsma, R Hiemstra, J Kreeft, A Palha, P Rebelo, and D Toshniwal. The 

geometric basis of mimetic spectral approximations. Proceedings oflCOSA- 

HOM 2012-2013 (this issue), 2012. 
[7] Marc Gerritsma. Edge functions for spectral element methods. In Spectral 

and High Order Methods for Partial Differential Equations, pages 199-207. 

Springer, 2011. 
[8] E Kanso, M Arroyo, Y Tong, A. Yavari, Marsden J.E., and M. Desbrun. On 

the geometric character of stress in continuum mechanics. Mathematik und 

Physik, 2007. 
[9] Jasper Kreeft and Marc Gerritsma. Mixed mimetic spectral element method for 

stokes flow: A pointwise divergence-free solution. Journal of Computational 

Physics, 240:284-309, 2013. 
[10] Jasper Kreeft, Artur Palha, and Marc Gerritsma. Mimetic framework on 

curvilinear quadrilaterals of arbitrary order. Arxiv preprint arXiv: 11 11.4304, 



page 69, November 201 1. URL http : //arxiv. org/abs/1111 . 4304 



[11] A PaUia, P Rebelo, and M Gerritsma. Mimetic Spectral Element solution for 
conservative advection. Proceedings of ICOSAHOM 2012-2013 (this issue), 
2012. 

[12] E. Tonti and Gruppo nazionale per la fisica matematica. On the formal struc- 
ture of physical theories. Istituto de matematica, Politecnico, 1975. 

[13] A Yavari. On geometric discretization of elasticity. Journal of Mathematical 
Physics, 2008. 



